Cluster Algorithm for hard spheres and related systems 
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In this paper, we present a cluster algorithm for the simulation of hard 
spheres and related systems. In this algorithm, a copy of the configuration 
is rotated with respect to a randomly chosen pivot point. The two systems 
are then superposed, and clusters of overlapping spheres in the joint system 
are isolated. Each of these clusters can be "flipped" independently, a process 
which generates non-local moves in the original configuration. A generaliza- 
tion of this algorithm (which works perfectly well at small density) can be 
successfully made to work at densities around the solid-liquid transition point 
in the two-dimensional hard-sphere system. 
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Since the 1987 paper by Swendsen and Wang |l|] (and the subsequent paper of Wolff 
the simulation of Ising or XY-tjpe systems close to the critical point has been much simpli- 
fied: just as the physical system, conventional Monte Carlo algorithms (for a review cf. [Q) 
suffer from critical slowing down, but the new algorithms overcome this problem and allow 
the calculation of thermodynamic quantities with great ease. 

One of the long-standing problems in classical statistical physics is the hard-sphere liq- 
uid [0. In two dimensions, the transition between the liquid and the solid order in the 
hard-sphere liquid has been the subject of unabating interest 0. There are several com- 
peting theoretical scenarios for the transition, and Monte Carlo work has been going on for 
more than 30 years P (for a review cf. 0). However, the conventional local- move Monte 
Carlo simulations are greatly affected by the slowing down of the simulation around the 
transition. At present, the maximum size of the simulation box, which can be unequivocally 
thermalized, contains only of the order of 1000 particles ||]. Much larger simulations have 
been undertaken ^,TU\ and sophisticated data analysis has been performed [TD|. However, 



due to the fact that the probability distribution has not yet converged to its equilibrium 
value, these simulations are biased in a way which is very difficult to assess. 

In this paper, we present a cluster algorithm, which is applicable to the hard sphere sys- 
tem in any dimension, and which is easily generalized to incorporate an additional potential. 
The main idea of the algorithm is to rotate a copy of the "current" configuration, and to 
superpose this rotated copy with the original simulation box. Clusters are then isolated 
in the joint system. Each of the clusters is then fiipped independently, i. e. the spheres 
belonging to a cluster are moved from the rotated copy into the original configuration and 
vice versa. For concreteness, consider the fig. 1, which illustrates the algorithm: the original 
configuration Ci (fig. la) is made up of N spheres of radius r in a box of size (L^,, Ly) (in the 
figure = 9). In addition to the original configuration, we consider also a configuration C2 
(displayed in fig. lb), which is obtained from ci by a tt- rotation: We generate C2 by picking 
an arbitrary "pivot" point p = {px,Py) with < p^ < Lx,0 < Py < Ly (In the example, 
Px = 0.52 Lx,Py = 0.53 Ly). We then rotate ci around p of an angle vr to obtain C2 [0. The 
choice of the pivot is the essential Monte Carlo element of the algorithm. 

The two configurations Ci and C2 are then superposed as shown in fig. Ic, where they 
form clusters of overlapping spheres 0]. Two types of clusters are possible: "even" clusters, 
made up of an equal number of spheres in ci and C2, and "odd" clusters, in which the numbers 
differ. In fig. Ic, e. g., the cluster "II" is even, while cluster "I" is odd. 

Generally, clusters appear in pairs (such as I and IV), with the possible exception of a 
single "even" cluster which is symmetric around the pivot (III). 

It is now easily seen that we may "fiip" the clusters, i. e. interchange between ci and C2 
the spheres belonging to a cluster. We are interested in performing a canonical simulation 
in which the individual numbers of spheres both of ci and C2 have to remain unchanged. 
We therefore choose to perform such fiips for individual even clusters, or for pairs of odd 
clusters, such as I and IV. The result of one such cluster flip (of clusters I and IV) is shown 
in fig. Id. Finally, we restrict our attention back to the updated configuration, c[ in the 
original simulation box. By inspection of fig. le, which shows €[, we see that picking p and 
fiipping clusters I and IV has achieved a nonlocal Monte Carlo move: spheres 8 and 9 were 
moved from the lower left corner to the upper right one, and the sphere 6 from the upper 
right corner to the lower left one. 
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It should be evident that - given an arbitrary pivot point - the flip satisfies the detailed 
balance condition, and constitutes a viable Monte Carlo move. To see this one simply needs 
to consider the "reverse" move (from fig. le back to fig. la), which has exactly the same 
probability to occur as the original one. 

Applying the same argument as above to an even cluster (such as V in fig. Ic), we 
notice that small even clusters generate only local moves. Due to the limited benefits of 
worrying about (small) even clusters, we usually exclude them from our considerations. 
Many generalizations are possible: It is evident that the spheres can have different radii etc] 
a potential can be taken into account in usual way, by calculating the Boltzmann weights 
of the proposed flip and the reverse one; furthermore, the angle of the rotation around the 
pivot can be chosen at will. This only introduces some effects far from the pivot, which can 
be eliminated. 

The simple algorithm which has just been described works perfectly well. At small 
density, the combined system of Ci and C2 breaks up into a large number of small clusters, 
which can be flipped independently. At higher density {i. e. above the percolation threshold 
of this combined system), there is a single percolating cluster (which it is useless to flip), and 
an algebraically decaying distribution of small clusters [0]. Just as in the Swendsen-Wang 
algorithm, there is a "magical" point, the percolation threshold (which, for the Ising model, 
corresponds to the Curie temperature At this point, the behavior of the algorithm 



is optimal. In our case, the "magical point" is the percolation threshold of the system of 
superposed conflgurations, which unfortunately lacks physical interest. In two dimensions, 
we flnd this percolation threshold to be situated at a density of p ~ 0.62, deflnitely lower 
than the densities in which we are interested (as usual 0, the density is deflned as the ratio of 
the number of spheres and the volume of the simulation box, p = N/V, normalized to 2/-\/3 
for the most compact state; in these units, the transition takes place around p ~ 0.9 P,p!0[1). 
Around the percolation threshold, the algorithm decorrelates the whole system by flipping 
a few large cluster, as in the Wolff algorithm . 

Close to the liquid-solid transition in the two-dimensional system, i. e. much above the 
percolation point, it is particularly difficult to flnd a sufficient number of small odd clusters, 
which generate the non-local moves. It is difficult to find odd clusters, but one rather often 
encounters configurations which almost constitute odd clusters, as the ones presented in fig. 
2, which are kept from fiipping by a few weak "links". We now present a stochastically 
correct trick which has allowed up to break up a large number of these weak links. 

For fixed but arbitrary e, we define an e— cluster as a set of spheres which may have an 
arbitrary number of e— links, i. e. links between a disk of the set and a disk of the boundary 
(in a different box), larger than 2 x r — e (so that the overlap between the spheres is smaller 
than e). In addition to having e— links, the e— cluster itself is held together by links which 
are not e— links, i. e. which are shorter than 2 x r — e. 

After isolating an e— cluster, we "freeze" the boundary, and perform a certain fixed 
number of nioc local Monte Carlo moves exclusively of the spheres in the e— cluster (without 
destroying it). Each time the number of e— links falls to zero, we have obtained a true cluster, 
which we fiip. It is easily shown that this exotic "dynamic of e— clusters" satisfies detailed 
balance, since we make sure that at each step both the initial and the final configuration 
are e— clusters. 

We have programmed the complete algorithm sketched in this paper. For small systems 



3 



(up to 14 spheres), we performed extremely long runs both of a standard MC algorithm 
and of the present one. We find identical probability distributions e. g. for the orienta- 
tional order parameter to a precision of 0.1%. There is thus very little room for doubt 



about the correctness of the present algorithm, and for programming errors in our actual 
implementation. 

It is easily understood that, at high density, the main workload of the algorithm consists 
in the determination of the percolating cluster. Since we never actually "flip" this cluster, 
most of the effort is thus spent in finding out what one does not want to do - a frustrating 
way of using CPU time . Only after discarding the percolating cluster, we have a chance of 
finding small odd clusters. A moment's thought suffices to understand that there is a faster 
way to find the small odd clusters. Consider the cluster IV the fig. Ic): it is evident that 
the nonlocal move can be performed only under the condition that, locally, it is possible 
to replace the sphere 6 by two other spheres. Whether there is at all enough space in a 
given neighborhood to replace n spheres by n + 1 can (for small n = 0, 1, 2) be decided by 
a local analysis which, in addition, is independent of the pivot (this remark is pertinent to 
clusters and e— clusters). An improved algorithm of the kind presented in this paper thus 
first isolates the loci at which n spheres may be replaced by + 1 (for small n = 0, 1,2). 
Once this analysis is done, one chooses randomly pivots, and then does the - now trivial - 
cluster search, in regions which have survived the screening stage. Since the local analysis 
just described can be reused a large number of times (and updated, once we have flipped 
a cluster, or e— cluster), the cluster search is much simplified. The final algorithm is thus 
quite efficient in generating non-local moves. 

Let us finally stress that the method presented in this paper is very general, and may 
be applied to a large number of systems. The specific application to the hard-sphere liquid 
stands out as prototypical, and we hope that the algorithm will be helpful in elucidating the 
order of the transition, in the measure of correlation functions etc. Work along these lines 
is in progress. 
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Figure Captions 



1. The algorithm presented in this paper performs nonlocal moves ( a) — > e) by consid- 
ering a randomly rotated copy ( b) ) of the actual configuration (a)), a) and b) are 
superposed ( c) ) and clusters are isolated ( c)) and flipped ( d)) in the superposed 
configuration. 

2. We present a trick which allows to flip not only clusters (as in fig. |l|), but also 
e— clusters. 
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